############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
base::library(stargazer)
base::library(DataCombine)
base::library(haven)
base::library(nnet)
base::library(ggeffects)

load(here("Data", "IndLevel_Data.RData"))
glimpse(ces_data)

table(ces_data$outcome_choice)
ces_data <- ces_data %>% 
  mutate(fac_outcome = factor(outcome_choice,
                              levels = c("Rep", "Dem", "Abst"),
                              labels = c("Rep", "Dem", "Abst"))) %>% 
  mutate(rural_urban = as.factor(rural_urban))


ces_data <- ces_data %>%
  mutate(party_id3 = ifelse(pid7 == 1 | pid7 == 2 | pid7 == 3, 
                            "Democrat", NA),
         party_id3 = ifelse(pid7 == 4, "Independent", party_id3),
         party_id3 = ifelse(pid7 == 5 | pid7 == 6 | pid7 == 7, 
                            "Republican", party_id3)) %>% 
  mutate(party_id3 = factor(party_id3,
                            levels = c("Democrat", "Independent",
                                       "Republican")))

#log_death, death_rate, 
summary(ind_model <- multinom(fac_outcome ~ death_rate*party_id3 + unemp_rate + defl_pcincome +
                                log_pop + rural_urban + defl_pcincome + #pid7 +
                                gender + black + hispanic + asian + other_race + 
                                log_age + college + income + year_elec, 
                              data = ces_data, weights = weight_cumulative))
#summary(ind_model)
#library(sandwich)
#coeffs_cl <- coeftest(ind_model, vcov = vcovCL, cluster = ~ fips)
#coeffs_cl

summary(ces_data$death_rate)
model_pid3 <- ind_model %>% 
  ggeffect(terms = c("death_rate[0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120]", 
                     "party_id3"), 
           ci.lvl = 0.90)

glimpse(model_pid3)
table(model_pid3$response.level)

model_drugs <- data.frame(model_pid3)
glimpse(model_drugs)
#rep_vote <- model_drugs %>% filter(response.level == "Rep")


vote_rep <- model_drugs %>% 
  filter(response.level == "Rep") %>% 
  mutate(group = factor(group,
                           levels = c("Democrat", 
                                      "Independent",
                                      "Republican"),
                           labels = c("Democrats", 
                                      "Independents",
                                      "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability",
       title = NULL) +
  ggtitle("A. Republican Vote") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), lim = c(0, 1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))

vote_rep

vote_dem <- model_drugs %>% 
  filter(response.level == "Dem") %>% 
  mutate(group = factor(group,
                        levels = c("Democrat", 
                                   "Independent",
                                   "Republican"),
                        labels = c("Democrats", 
                                   "Independents",
                                   "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability",
       title = NULL) +
  ggtitle("B. Democratic Vote") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), lim = c(0, 1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "none",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))

vote_dem

abst <- model_drugs %>% 
  filter(response.level == "Abst") %>% 
  mutate(group = factor(group,
                        levels = c("Democrat", 
                                   "Independent",
                                   "Republican"),
                        labels = c("Democrats", 
                                   "Independents",
                                   "Republicans"))) %>% 
  ggplot(aes(shape = group, color = group)) + 
  geom_pointrange(aes(x = x, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .85#, position = position_dodge(width = .5),
                  #color = "darkslategray"
  ) + 
  theme_light() + 
  labs(x = "Overdose Death Rate", y = "Predicted Probability",
       title = NULL) +
  ggtitle("C. Abstention") +
  scale_shape_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c(15, 16, 17)) +
  scale_color_manual(name = "Party Identification:",
                     labels = c("Democrats", "Independents", "Republicans"),
                     values = c("blue", "gray40", "red")) +
  scale_x_continuous(breaks = c(0, 20, 40, 60, 80, 100, 120), lim = c(-1, 121)) +
  #scale_y_continuous(breaks = c(0, 0.1, 0.2, 0.3), lim = c(0, 0.27)) +
  scale_y_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), lim = c(0, 1)) +
  theme(axis.text.x = element_text(hjust = 0.5, size = 11),
        axis.text.y = element_text(size = 11),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 12),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))

abst

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(vote_rep, abst, vote_dem, cols = 2)

